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I. INTRODUCTION 

Helium is the second element in the periodic table, 
as well as the second most abundant chemical element 
in the universe. It is a major constituent of both stars 
and giant planets, and it is involved in several different 
nuclear fusion reactions. Its high-pressure physical prop¬ 
erties are therefore rather important for many different 
branches of natural science. The low-temperature be¬ 
haviour of helium is well-known but rather peculiar: it 
stays liquid up to the absolute zero while becoming super¬ 
fluid, quantum-freezes under pressure, and demonstrates 
quantum-crystal behavior in the solid phase. These 
quantum effects are a consequence of a relatively small 
mass of the He atom (« 7296 electron masses for 4 He) 
and weak interatomic interaction. Solid helium is a close- 
packed atomic crystal with nearly spherically symmetric 
He atoms. 4 He has hexagonal close-packed (hep) struc¬ 
ture everywhere except for small body-centered cubic 
(bee) and face-centered cubic (fee) regions near the melt¬ 
ing line on the (p, T) phase diagram!. The c/a ratio of the 
hep structure is close to the ideal value 8/3 « 1.633^. 


While helium has been an object of intensive experi¬ 
mental study for more than a century, the high-pressure 
experimental breakthrough happened in the last few 
decades due to the invention of diamond anvil cells. 
The elastic moduli of hep helium up to 32 GPa, as 
well as sound velocities, Poisson’s ratio (PR) and elastic 
anisotropy parameters, have been measured experimen¬ 
tally by Zha et alA The results were somewhat unex¬ 
pected. First, a significant anisotropy of elastic proper¬ 
ties was found. The problem of the elastic anisotropy 
of helium is important for the high-pressure experimen¬ 
tal techniques, as helium is frequently used as a quasi¬ 
hydrostatic medium!. Second, the Poisson’s ratio was 
found to decrease with increasing pressure, which is a 
rather unusual behavior, as for most solids PR increases 
with pressure, approaching 1/2 at megabar pressures. 
Similar decrease of PR with pressure has been observed 
in solid hydrogen!. Such anomalous PR behavior of He 
was often thought to be a quantum zero-point vibration 
(ZPV) effect, however, as we show in the present pa¬ 
per and in Ref. @, this is not the case. A theoretical 
calculation of elastic moduli by Nabi et al£ using den- 
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sity functional theory (DFT) in local Airy gas (LAG), lo¬ 
cal density approximation (LDA) and generalized gradi¬ 
ent approximation (GGA) soon followed the experiment. 
The calculated elastic moduli were in a reasonably good 
agreement with experiment. Unfortunately, the authors 
of Ref. H did not calculate Poisson’s ratio, and did not 
study the elastic anisotropy in any detail either. 

The goal of the present paper is to clarify the two is¬ 
sues introduced above. We calculate the elastic moduli 
of hep He as a function of pressure using two complemen¬ 
tary methods: DFT-GGA and semi-empirical (SE) po¬ 
tentials, specifically focusing on the Poisson’s ratio and 
elastic anisotropy parameters. We extend our GGA cal¬ 
culations into the metallic phase (up to p = 30 TPa) in 
order to check the effect of terapascal pressures and the 
metallization transition on the elastic properties of he¬ 
lium. We also study the effect of zero-point vibrations 
in the Debye approximation on the physical quantities in 
question (for 4 He) in order to determine whether quan¬ 
tum effects play any significant role at high (p > 10 GPa) 
pressures. 

The paper is organized as follows. In section [TT] we 
outline the two approached for calculating total energy 
(DFT-GGA and SE), and then give a brief introduction 
to the elasticity of hexagonal crystals, the algorithm of 
calculating elastic moduli of an hep crystal^, and define 
physical quantities used in the present paper. In section 
mu we present results of our calculations. 


II. METHOD 
A. Total energy calculations 

Helium consists of electrons and nuclei, and due to 
the relatively small mass of the latter, their quantum 
zero-point motion cannot be ignored in general. Many 
different numerical methods have been applied to solid 
helium^ - —. The electronic subsystem can either be de¬ 
scribed from first principles, like in density functional 
theory (DFT)-based methods, or replaced by empirical 
pairwise or n-body interactions between nuclei. The 
quantum and thermal motion of the nuclei can be ana¬ 
lyzed either in harmonic approximation, or using anhar- 
monic approaches, such as diffusion Monte Carlo (DMC) 
and other Monte Carlo methods^. 

For the elasticity calculations pair potentials are not 
sufficient. We have opted to use DFT in the generalized 
gradient approximation (GGA) as our main method. We 
have also used pair and 3-body empirical potentials for 
the p < 100 GPa range. For the quantum zero-point 
vibrations we have used a simple harmonic Debye ap¬ 
proximation (see below). While such approach is rather 
crude, and the use of harmonic approximations for He 
atoms has been recently criticized in Ref. E we chose it 
as it is computationally cheap and consistent with a full 
DFT treatment of the electronic subsystem. 

For the density functional theory calculations we have 


used the all-electron full-potential linear muffin tin or¬ 
bital (FP-LMTO) code RSPti^ with the GGA functional 
of Perdew, Burke, Ernzerhof (PBE)i^. The basis set in¬ 
cluded Is, 2p and 3d electrons of helium, with two LMTO 
basis functions with kinetic energies —0.1 Ha and +0.1 
Ha respectively per each atomic orbital. 847 ^-points in 
the Brillouin zone have been used. 

For the semi-empirical calculations we have used the 
pair and triple SE potentials described in Ref. UH. They 
include Aziz pair potential in the Silver a-Goldman form 
and the three-body potential in the Slater-Kirkwood 
form. These are exactly the potentials used in our pre¬ 
vious works^^ d 6 - 17 on helium. The cutoff radii R 2 = 
50.2a and R 3 = 10.2 a were used for pair and triple forces 
respectively, with a being the lattice constant. 

Both DFT and SE calculations were performed for 
zero temperature, and the zero-point vibrations were ne¬ 
glected at first. The effect of ZPV was later accounted 
for in the Debye approximation. In contrast to Ref. H, 
we found the pressure-dependent equilibrium c/a ratios 
as described in Refs. @ and H and used them for all our 
calculations. All our results are well converged with re¬ 
spect to the number of ^-points (GGA) and cutoff radii 
(SE) respectively. 

DFT-GGA and SE can be seen as complementary 
methods. Semi-empirical potentials usually work very 
well for low pressures, but fail for higher pressures where 
4-body and higher order n-body forces become impor¬ 
tant. For helium the threshold pressure is of the order 
of 100 GPa£^. GGA, on the other hand, can be inaccu¬ 
rate at low pressures due to poor description of the van 
der Walls (vdW) forces. It has been shownii that the 
effect of the vdW forces is rather small in the gigapascal 
pressure range. 

B. Elasticity under initial pressure 

This subsection and the remaining part of section El 
deal with the elasticity theory for a hexagonal crystal un¬ 
der pressure and the method of calculating elastic mod¬ 
uli numerically. It includes the main results of Refs. 
SE and 19 , and other works. This material has never 
before been gathered in one place, therefore we decided 
to give a brief introduction to the method as a whole, 
presenting all relevant formulas and stressing some im¬ 
portant points. 

If a strain is applied to an elastic medium, each point 
r of the medium is shifted to a new position r'(r), and 
we can define tensor uij as 


where e^- = (u^ + Uji )/2 and cOij = (uij — Uji )/2 
are the symmetric strain tensor and the antisymmet¬ 
ric rotation tensor respectively, and the tensor indices 
= 1,2,3 number three Carthesian (not crystal) co¬ 
ordinates. Summation over repeated tensor indices is 
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assumed. In the present paper we do not consider ro¬ 
tations, so we always assume that = uji = eij and 
ojij = 0. Note that we are not using the Lagrangian 
strain tensor rjij = (u^ + Uji + UkiUkj)/2 in the present 
paper. The difference between and rjij is important 
under external pressure. 

The volume of the strained medium is 

V' = V det (Sij + Uij ) = V + AV + 0(u 3 ), (2) 

where 


uniform strain (i.e. the lattice geometry is not allowed 
to change as the volume changes): 

(8) 

V'=V, fixed geometry 

It can be shown that Ky (p) = Cukk (p) /9 and it does not 
depend on p explicitly. Reuss bulk modulus is defined 
as the bulk modulus under the uniform stress (i.e. the 
lattice geometry is changed as the pressure changes): 


K v 


V 


d 2 E 


dV' 2 


AV = V 


l, n2 1 

+ ofaii) - o 1 


( 3 ) 


is the change of volume up to the second order in , and 
we have used the identity det A = exp Tr log A to expand 
V' in powers of uij. 

The elasticity theory for a medium under external 
stress erf? is rather non-trivia l 18 i 2Q and there is no 
straightforward generalization of the zero-pressure stiff¬ 
ness tensor Cijki • However, theory simplifies for the case 

of the isotropic external pressure = —pSij. The 
stress-strain relation up to the first order in isi£ 

<?ij = Cijki ( p)uki ~ pSij , (4) 


K r = -V 


d?E 

~dV ?2 


= V 

V'=V, relaxed geom. 



( 9 ) 

Kr is equal to 1/Sukk(p), where Sijkiip ) is the compli¬ 
ance tensor, the inverse of Cijki (p ), and again it does 
not depend on p explicitly. The shear moduli are calcu¬ 
lated from Eq. © using strains which are isochoric 
(volume-conserving) exactly, or in the second order of 
at least, so that AV = 0 in Eq. © and the elastic 
energy is proportional to the shear modulus in question 
times u 2 , without any additional terms proportional to 
pu 2 . 


Hexagonal close-packed (hep) crystal lattice has lattice 
vectors 


where the rank-four pressure-dependent stiffness tensor 
Cijki (v) (tensor of elastic moduli, called c a p crr in Ref. 
[Ts) has the same symmetry as the zero-pressure Cijki , 
namely 

Cijki = Cjiki — Cijik = Ckuj • ( 5 ) 


" 1" 


'-1/2' 


"o" 

0 

,a 2 = a 

V3/2 

,a 3 = c 

0 

0 


0 
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where a and c are two pressure-dependent lattice con¬ 
stants. Atom positions in the undeformed hep lattice are 


The elastic energy density up to the second order in u\ 
^18 


e =y = 2 C ijkl(p)uij u kl ~ P—, (6) 

where AV is defined in Eq. ([3j) . The terms of the order 
u 2 in AV are important, as they are of the same order 
as the first term in Eq. <©• Note that the energy density 
e is defined with respect to the undeformed volume V, 
not V'. Total energy E gives adiabatic stiffness tensor, 
while for the isothermic one the Helmholtz free energy 
F = E — TS should be used instead. In this paper we 
limit ourselves to the case T = 0, so there is no distinction 
between the two. The equation of motion of the elastic 
medium up to the first order in i 


d 2 Ui 

p dt 2 


= Cijki (p) 


djuk 

dxjdxi ’ 


( 7 ) 


where iq = x\ — Xi. 

When calculating elastic moduli numerically from Eq. 
([6]) one must be careful with the —pAV/V term, as it 
includes u 2 terms. One way of addressing the problem is 
to calculate bulk and shear moduli separatelySiSl. Voigt 
bulk modulus is defined as the bulk modulus under the 


(0,0,0), (y,i) an 

in crystal coordinates. When a uniform strain is 
applied to the medium, any point Xi changes to x\ = 
(Sij + Uij)xj , i.e. atoms form a deformed crystal lattice 
with new lattice vectors 

a i = (fiij 4 " u ij) a j • ( 12 ) 

Special care must be taken when applying strains to a 
lattice with more than one atom per unit cell. The ex¬ 
pression m determines only the change of three lattice 
vectors under strain, but tells nothing about positions 
of atoms within the unit cell. Such positions are not 
determined by macroscopic elasticity theory and must 
be allowed to relax in total energy calculations. If they 
change linearly in under strain, this affects the calcu¬ 
lated elastic constants. In other words, to obtain correct 
result we must minimize the total energy with respect to 
all atomic positions for each applied finite strain. Hep 
lattice has two atoms per unit cell, and their positions 
are fixed by symmetry for the undeformed lattice, how¬ 
ever for certain strains (like the orthorhombic strain, see 
below) they indeed change linearly in u^. 
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In Voigt notation 


The corresponding strain (found from Eq. fl3J)) is 


(11, 22,33, 23,31,12) (1, 2,3,4, 5, 6 ) (13) 

the stiffness tensor of a hexagonal crystal is specified by 
five independent elastic constants Cn, Ci 2 ,Ci 3 , C 33 and 

C 44 : 


Cij 


Cil C13 0 0 

C 12 Cn Cis 0 0 

Ci3 C13 C33 0 0 

0 0 0 C44 0 

0 0 0 0 C44 

0 0 0 0 0 


0 
0 
0 
0 
0 

>(c n — C 12 ) _ 


• ( 14 ) 


C. Calculation of elastic constants 

If we have any method of calculating the energy of the 
crystal for given volume V, c/a ratio and strain Uij, like 
DFT or SE, we can calculate the five elastic constants 
numerically by applying several independent strains in 
Eq. ([ 6 |) . An efficient way to do this for the hep lattice has 
been proposed by Steinle-Neumann et al£. As explained 
above, bulk and shear moduli are calculated separately. 

First we calculate the total energy E(V) as a function 
of volume (in practice the molar volume V mo i = VNa/N 
is used, where N is the number of atoms, and Na is Avo- 
gadro’s number). For each volume the energy minimum 
with respect to the ratio c/a is found, and E(V) is de¬ 
fined as the energy of this minimum. The ideal c/a ratio 
(corresponding to the close-packing of hard spheres) is 
y/8/3 « 1.633, and for helium the c/a ratio is rather 
close to the ideal valued. The equation of state (EOS) 
p(V) and the Reuss bulk modulus Kr are found as 

_ dE _ dp _ d?E 

P dV' R V dV dV 2 ' ( 1 

In practice, the energy E(V) is approximated (via least 
square fit) by the Rose-Vinet equation of stated, which 
allows for an accurate numerical differentiation. For 
GGA, we were unable to find a single Rose-Vinet fit 
which would be accurate in both GPa and TPa pressure 
ranges, therefore we had to use two different parametriza- 
tions. 

By using a diagonal stress (Jij = — (p + A p)Sij, which 
is an infinitesimal change of pressure, we can show that 
the Reuss bulk modulus is 


S- Q 

Kr c? 


(16) 


A p 

0" 


C 33 — Cl 3 
0 
0 


0 

C 33 — Ci 3 

0 


Cn + C 12 — 2 C 13 
(19) 

Eq. (fl6l) is the first equation we use for determining 
five elastic constants. From Eq. m we can find the 
logarithmic derivative of c/a with respect to volume 


R = _dlog(c/a) 


= ^(C 33 + C 13 -C n -C 12 ). ( 20 ) 


d log]/ C s 

This is the second equation we need. An analytic 
parametrization of c/a(y mo i ) is used for numerical dif¬ 
ferentiation as usual. 

From now on we are going to use only exactly isochoric 
(volume-conserving) strains as explained above. In order 
to calculate Cs, Ref. M used an isochoric c/a-changing 
strain, however it is not necessary, as we can use an equiv¬ 
alent formula 


Cs 


9 /c \ 2 d 2 E(V, c/a) 


-(-X 

2 V \aJ 


d ( c/a) z 


( 21 ) 


taken at the equilibrium c/a. We use the energies 
Efy^c/a) calculated previously when we looked for the 
equilibrium c/a for each volume and approximate it by a 
fourth order polynomial of c/a. Equation (fl 8 |) with Cs 
from Eq. m is the third equation for the five elastic 
moduli. We now have three equations for three variables 
Cn + C 12 (in this combination only), C 13 and C 33 . They 
can be solved to obtain 

C u +C 12 = 2 K r + ^-(2 R-1) 2 , ( 22 ) 


C 1s = K r +^(2 R-1)(R+1), (23) 


C 3 3 = K r + ^(R+ l) 2 . 


(24) 


We need two more isochoric strains to find C 44 and 
Cqq = (Cn — Ci 2 )/ 2 . C 44 is found from the monoclinic 
strain 


0 


0 ^ 0 
t 0 0 


(25) 


where 


Q = C 33 (Cn + C 12 ) — 2 Ci 3 , 


(17) 


with e = 2 C 44 C + O(C). Alternatively, C 44 for the hep 
lattice can be obtained from the calculated Raman 
frequency using the formula^ 


(18) 


C 44 = 


m 


4a/3 cl- 




29 ’ 


Cs = Cn + C12 + 2C33 — 4 Ci 3 . 


(26) 
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as we did previously in Ref. Uzl The resulting values for 
C 44 , obtained using these two approaches agree within 
1%. For Cqq we use the orthorhombic strain 


"tO 0 ' 

0 -t 0 

. 00 ^.’ 


(27) 


which gives e = 2C^t 2 ± 0(t 4 ) = (Cn — Cy^t 2 ± 
0(t 4 ). Note that for this strain the atomic positions 
change under strain. The position of atom 2, which 
is (2/3,1/3,1/2) for the undeformed lattice, becomes 
((l±A)/ 2 ,A,l/ 2 ), with parameter A depending linearly 
on t as A = 1/3 — gt ± 0(t 2 ). It means that we have 
to relax the parameter A for each finite strain t in or¬ 
der to calculate the correct elastic energy. In our cal¬ 
culations of C 44 and Cqq we have used five t-points 
t = 0 , ±0.005, ± 0.010 (for very high pressures smaller val¬ 
ues of t were used) and interpolated the energy E(t) with 
a fourth order polynomial in order to find the coefficient 
before t 2 . The results obtained this way are virtually 
independent on the particular values of t used. Such ap¬ 
proach is vital for ± 66 , while for C 44 the E(t) dependence 
is almost exactly quadratic in t in a wide range of t. With 
C 44 and ±66 calculated we can finally determine all five 
elastic constants of the hep lattice. The quantities C 44 , 
Cqq and Cs/6 are three different shear moduli for three 
independent pure shear (isochoric) deformations. They 
are all equal for an isotropic solid. 


D. Anisotropy parameters 


The three acoustic sound velocities (one compressional 
and two shear ones) of the hexagonal lattice are &21 


P v s 1 


2 A -\- B 
P V P = o > 


2 _ Cll Cl2 .2 n I n 


sin 0 ± C 44 cos #, 

A-B 


where 


P v s 2 = 


A = Cn sin 2 0 ± C 33 cos 2 0 ± C 44 , 


(28) 

(29) 

(30) 

(31) 


B 2 = [(Cn — C 44 ) sin 2 9 ± (C 44 — C 33 ) cos 2 0] (32) 


±(Ci 3 ± C 44) 2 sin 2 ( 2 $), 

and 0 is the angle between wave vector q and z. 

The elastic anisotropy of a hexagonal crystal can be 
described by the anisotropy parameters of these three 
acoustic waves^ 


C 33 

C 


ASi = 


Cn ± C 33 — 2 C 13 

4C 44 


Cs ± 2 C 66 

8C44 


(34) 


AS 2 


2C 44 

Cn — C12 


C 44 

Cqq 


(35) 


For an isotropic medium this quantities are equal to 
one. Note that for a cubic crystal A P is equal to 
unity, and A Si is the single anisotropy parameter (Cn — 
C12)/(2C44) of the cubic crystal. AS 2 is ambiguous, as 
there is no condition Cqq = (Cn — Ci2)/2 for the cubic 
symmetry, but instead C44 = Cqq. 


E. Aggregate properties: Voigt and Reuss 
approaches 


Aggregate description replaces an actual crystal with 
an effective isotropic elastic medium, described by an av¬ 
erage bulk modulus K and an average shear modulus G, 
with 

k ~I g ) s 


ij&ki + C (SikSji ± SuSjk) • (36) 



It is a good approximation for polycrystalline solids and 
mixtures. The compressional and shear sound velocities 
are 





(37) 


This is often written as v 2 P = v 2 B ± with vb = {K/p) 1 
being the bulk (hydrodynamic) sound velocity. It has no 
direct physical meaning for a crystal, however it corre¬ 
sponds to the sound velocity of the liquid phase, as can 
be seen at the liquid-solid transition in Ref. @. The Pois¬ 
son’s ratio is defined as: 


13K-2G 1 3v% - 2vl 1 v 2 P - 2vl , N 

aS 2jKTG =2^4=2-^- <38) 

In order to find the values of K and G consider a poly¬ 
crystalline solid consisting of grains of all possible orien¬ 
tation, and use Voigt and Reuss estimates of its elastic 
moduli. In the Voigt approach a uniform strain field Uij 
(of either uniform compression or a shear type) is applied 
to the mixture, and K and G are found from the averaged 
elastic energy density. For a polycrystal this means aver¬ 
aging the elastic energy over all possible orientations of 
the crystal relative to the strain field . Reuss approach 
uses uniform stress instead. The bulk moduli defined in 
this way are equal to Ky and Kr defined above. For a 
hexagonal lattice the bulk and shear moduli arei^ 


K v = I (2C n + C 33 + 2Ci 2 + 4 C 13 ), I< R = ± (39) 

9 C s 


Gy 


30 


(Cs ± I2C44 ± 12 C 66 ) 5 


(33) 


(40) 
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5 QC^Cqq 

2 SKvC^Cqq + Q(C44 + Cqq) 


(41) 


Finally we use Voigt-Reuss-Hill average scheme to obtain 
K and G : 


K = 


Kv+K r 

2 


G = 


Gy + Gr 
2 


(42) 


F. Zero-point vibrations 


In our calculations we have treated zero-point vibra¬ 
tions within the framework of the Debye model. As dis¬ 
cussed above, this approach is by no means exact, espe¬ 
cially at low pressures, but it gives us a good estimate 
of the effect of ZPV on the elastic properties. Using K 
and vs obtained in the absence of ZPV, we calculate the 
Debye temperature 


T d = 


Jzb 


zero-point energy 


Vmol (± 2_ 

187 T 2 Na \Vp Vg 


E zp = -NksTD, 


-1/3 


(43) 


(44) 


and the corresponding contributions to pressure and the 
bulk modulus: 


A p = — 


dE zp 

~dV~ 


9 AT j dT D 

s NAkB dv^r 


(45) 


A T _ dAp 9 Tr _ d 2 Tr> 

AK = -Vmol-— - yVmolNAkB 2 ‘ (46) 

dVmol 8 dV 2 ol 

We find the volume derivatives of Tr> analytically, us¬ 
ing the Rose-Vinet expression for K and a simple 
parametrization for vs- Below we always compare results 
with and without ZPV in order to measure the impor¬ 
tance of quantum effects on each physical quantity. 4 He 
isotope was assumed for all our ZPV calculations. 


III. RESULTS AND DISCUSSION 

A. Equation of state, metallization and the c/a 
ratio 

The calculated equations of state (EOS) p(V mo i ) of the 
hep helium are presented in Fig. |T] and compared with 
the experimental data£. The four curves correspond to 
our four approaches: SE and GGA, with and without 
ZPV. EOS can be viewed as an auxiliary quantity in our 
calculations, as it is used to calculate Reuss bulk mod¬ 
ulus. It is also used to change variables from V mo i to p, 
which is done to obtain p-dependent quantities presented 
in all figures of the present paper. For this we always used 
the respective EOS for each of the four approaches, e.g. 



FIG. 1. (Color online) Equation of state of hep helium. Inset: 
high-pressure region. The vertical line indicates the GGA 
metallization point. 


a GGA+ZPV EOS was used for GGA+ZPV elastic mod¬ 
uli, but SE (No ZPV) EOS was used for SE (No ZPV) 
elastic moduli. 

The semi-empirical potentials underestimate the pres¬ 
sure p for a given V mo i, but zero-point vibrations im¬ 
prove the agreement with the experiment significantly. 
GGA overestimates the pressure, and the inclusion of 
ZPV makes things worse. The effect of ZPV is more pro¬ 
nounced in the SE approach, as the SE pressure for a 
given Vmol is smaller compared to GGA. 

The inset of Fig. [T| shows the high-pressure EOS 
(GGA with and without ZPV) for pressures up to the 
metallization point and above. The metallization takes 
place at V mo i = 0.228 cm 3 /mol in our GGA calculations^ 
(p = 17.48 TPa from our GGA+ZPV EOS, or p = 17.08 
TPa without ZPV). This point is shown as the verti¬ 
cal line in the inset of Fig. [lj Our metallization vol¬ 
ume and pressure are in good agreement with previ¬ 
ous GGA result s 12 ! 25 . Note, however, that GGA seri¬ 
ously overestimates the metallization volume (by about 
20% for He according to the diffusion Monte-Carlo and 
GW studie s 12 ! 25 ), thus underestimating the metalliza¬ 
tion pressure. Moreover, it has been recently showni^ 
that vibrational degrees of freedom further increase the 
metallization pressure, but these effect cannot be repro¬ 
duced in harmonic approximation. Ref. [H gives the 
value p = 32.9 TPa. Such questions are beyound the 
scope of the present paper. We use only GGA with ZPV 
in harmonic Debye approximation, which results in an 
underestimated metallization pressure of 17.48 TPa. 

Our calculated c/a ratios have been presented previ¬ 
ously in Refs. Q and 0. SE and GGA give somewhat 
different c/aiVmoi) curves, but both methods give lattice 
distortions S = cja — ^/8/3 of the order of 10 -3 at the 
pressures up to 150 GPa. 5 is negative for p > 13 GPa 
in both approaches. For the TPa pressures the magni- 
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FIG. 2. (Color online) Five elastic constants of hep helium 
versus pressure. 

tilde of the negative 5 obtained with GGA grows and it 
reaches a sharp minimum S ~ —0.05 at the metallization 
point Vmoi = 0.228 cm 3 /mol (Ref. d). 

B. Elastic moduli and the elastic anisotropy 

The five elastic moduli for the pressures up to 150 GPa 
are presented in Fig. [2j Again, we compare our four the¬ 
oretical approaches with the experimental data of Zha et 
alA. Note that all our calculations were done for T = 0, 
while the experimental elastic moduli were measured at 
room temperature T = 300 K, which can be one of the 
reasons for the theory-experiment discrepancies. Both 
SE and GGA are in reasonable agreement with the exper¬ 
iment, although neither method gives a perfect quantita¬ 
tive match. GGA seems to be more consistent of the two. 
Our GGA results without ZPV are close to the GGA and 
LAG results obtained in Ref. d- The DFT-GGA errors 
mainly stem from the inaccuracy of the GGA functional 
itself, we have checked that the parameters of the FP- 
LMTO calculations (number of Appoints, tail energies, 
etc.) have minimal effect on the results. It is unlikely 
that the errors of the SE approach (at least for the pres¬ 
sures p < 50 GPa) are caused solely by the neglect of the 
4-body and higher order n-body forces. The most likely 



FIG. 3. (Color online) Elastic moduli (GGA+ZPV only) of 
hep helium under terapascal pressures. Symbols are the calcu¬ 
lated data points. The vertical line indicates the GGA met¬ 
allization point. Cqq = {C ii — Ci 2)/2 and Cs/6 are also 
presented. 


reason for the SE-experiment differences (apart from the 
temperature effect) is that the elastic moduli are sensi¬ 
tive to the particular parametrizatiorJ^ of the pair and 
triple forces. The effect of ZPV, while noticeable, does 
not affect the elastic moduli in any drastic way. 

In Fig. [3] the calculated elastic moduli (GGA+ZPV) 
are presented for the terapascal pressure range. The two 
additional shear moduli Cqq and Cs /6 are also shown. 
Significant features are seen at the GGA metallization 
point p = 17.48 TPa. They are most likely finite discon¬ 
tinuities (steps), however, critical behavior at the metal¬ 
lization transition is beyond the scope of the present pa¬ 
per and would be difficult to investigate accurately with 
the methods we employ. Of the three elemental shear 
moduli only Cs/6 has a large jump (« —5%) at the met¬ 
allization point, while the behavior of C 44 and Cqq is 
relatively smooth. Remember that Cs is the measure 
of stiffness of the crystal with respect to the isochoric 
change of the c/a ratio, so it is not unexpected that the 
behavior of Cs is irregular at the point where c/a has a 
sharp minimum^. Note that Cs(p) also behaves rather 
nonlinearly in the metallic phase. These irregularities 
in Cs affect the five elastic constants, most notably C 13 
and C 33 , via Eqs. (f22l) - (|24|) . Hep lattice is dynamically 
stable for the whole pressure range considered, i.e. the 
conditions^ C44 > 0 , C\\ > IC 12 I, C 11 C 33 > (C 13) 2 and 
c^{c u + c 12 ) > 2 c ; 2 3 are fulfilled, which ensures that 
the single-crystal sound velocities (IZSJ-® are real. 

All our calculations were done for the relaxed c/a ra¬ 
tio, and the effect of the volume dependence of c/a has 
been taken into account through the parameter R. In 
order to test the importance of the c/a distortion for the 
elastic moduli we have also performed calculations with 
R = 0, like in Ref. [§|. Our results (not shown) indicate 
that the difference is barely visible for the pressures up 
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FIG. 4. (Color online) Elastic anisotropy parameters of hep 
helium versus pressure. 

to 150 GPa, however, it becomes significant (about 5%) 
for the TPa pressure range, and it noticeably affects the 
features near the metallization point. We took great care 
in choosing the parametrization of the function c/a(Vmoi) 
that reproduces the sharp minimum in c/a well. 

The elastic anisotropies A P, A Si and AS 2 are pre¬ 
sented in Fig. 0J The compressional anisotropy A P is 
close to the isotropic value of 1 in the experiment, how¬ 
ever both GGA and SE predict noticeable anisotropy 
A P > 1 . GGA gives a virtually pressure independent 
value A P ^ 1.14, while the semi-empirical A P grows 
from « 1.1 at low pressures to « 1.25 at 100 GPa. ASi 
is about 1.2 — 1.3 in the experiment, both GGA and SE 
overestimate it significantly, giving values of the order of 
1.6 — 1.7. The third parameter, AS 2 , is always less then 
one in both theory and experiment. SE and GGA, how¬ 
ever, give values closer to 1 than the experiment, thus 
underestimating the anisotropy. 

To summarize, GGA gives A P « 1.14, A Si « 1.7, 
AS 2 ~ 0.93 and the pressure dependences of these pa¬ 
rameters are rather small. SE approach, while agreeing 
well with GGA for low pressures, predicts much stronger 
pressure dependence of AP, ASi 5 2 , and different sign 
of dASi/dp. Ref. HI states that the three anisotropy 
parameters are nearly pressure-independent. While our 
GGA data fully confirms this result, our SE data behaves 
quite differently. In the absence of reliable experimen¬ 
tal data on the pressure dependences it is hard to say 



FIG. 5. (Color online) Elastic anisotropy parameters (up¬ 
per panel) and Cauchy violations (lower panel) of hep helium 
under terapascal pressures. Symbols are the calculated data 
points. The vertical line indicates the GGA metallization 
point. 


which behavior is more correct (see the discussion be¬ 
low, however). Both methods overestimate A P and ASi 
significantly, and also overestimate AS 2 (which under¬ 
estimates the anisotropy), thus neither approach agrees 
particularly well with the experiment. In particular, the 
relative errors in A P, A Si and AS 2 (computed relative 
to the experiment) are significantly larger than the errors 
in Cij. The reasons for such discrepancy are unknown. 
Since both GGA and SE agree with each other better 
than with the experiment, one can speculate that the 
difference in temperature between 0 K and 300 K might 
play at least some role. The large elastic anisotropy of 
He, while somewhat unexpected, is not in any way incom¬ 
patible with the high symmetry of the hep crystal and the 
nearly spherical shape of the atoms. Highly symmetric 
crystals have isotropic or nearly isotropic rank -2 tensor 
properties, such as conductivity, thermal expansion or 
dielectric permittivity. The stiffness, however, is a rank- 
4 tensor property, and it is well-known that even cubic 
crystals are not elastically isotropic. The large values of 
ASi are not very surprising, since this parameter is simi¬ 
lar to the single anisotropy parameter (C±i —C 12 ) /( 2 C 44 ) 
of the cubic crystal. 

The three anisotropy parameters at TPa pressures are 
presented in Fig. [5j upper panel. In the insulator phase 
anisotropy decreases with pressure, with three parame¬ 
ters becoming close to unity just before the metallization 
transition. This is an intuitively plausible behavior of 
atomic crystal becoming more isotropic under pressure. 
However, A P, A Si and AS 2 demonstrate strong fea¬ 
tures at the metallization point (especially ASi, which 
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FIG. 6 . (Color online) Cauchy violations of hep helium versus 
pressure. 


involves Cs), with anisotropies of two shear modes chang¬ 
ing sign at the metallization transition, and A P stay¬ 
ing close to one and showing nonlinear behavior in the 
metallic phase. One can understand the anisotropy pa¬ 
rameters in terms of the elastic moduli presented in Fig. 
[ 3 | For instance, let us analyze ASi. For an isotropic 
solid C 44 = Cq 6 = Gs/6. For helium at low pressures, 
however, C^/6 has almost twice the value of C 44 , giving 
a large AS\. At higher pressures Cs/6 is of the same 
order as C 44 and Cqq, and eventually becomes smaller 
than them in the metallic phase, which corresponds to 
A Si < 1. 

Strictly speaking, it is wrong to say that the anisotropy 
parameters are nearly pressure-independent. Indeed, 

the pressure-induced changes in Fig. [ 5 ] are rather dra¬ 
matic even if we consider the insulating phase only. The 
pressure-independence found in Fig. HJ and in Ref. 
is simply the result of the pressure scale of 150 GPa be¬ 
ing very small compared to the metallization pressure 
^17 TPa, which is presumably the only natural pres¬ 
sure scale in solid helium (if the quantum effects are dis¬ 
regarded). From this logic we can infer that GGA is 
more trustworthy than SE in determining pressure de¬ 
pendence of AP, A5i ? 2, since it gives reasonable results 
in TPa range of pressures, while the SE method breaks 
down at p ~ 100 GPa, which can act as a spurious pres¬ 
sure scale. This argument is far from infallible, however, 
as GGA is known to be unreliable for pressures p < 50 
GPa due to the poor description of the van der Walls 
forces. 

The Cauchy violations 3Ci2 — Cn—Ap and C13 — C 44 — 
2 p are presented in Fig. [6] They can be thought of as 
a measure of noncentral forces in a solid. Experimental 
3Ci2 — C11 — 4 p is reproduced well by SE and not so well 



FIG. 7. (Color online) Aggregate properties of hep helium 
versus pressure: bulk modulus K (upper panel), shear mod¬ 
ulus G (middle panel) and Debye temperature (lower panel). 
K and G are obtained using Voigt-Reuss-Hill averaging. 


by GGA, but the situation is reversed for C13 — C44 — 2 p. 
The Cauchy violations at terapascal pressures are shown 
in Fig. [ 5 j lower panel. Their behavior is mostly linear 
with mild kinks at the metallization point. 


C. Aggregate properties 

The Voigt-Reuss-Hill-averaged bulk and shear moduli 
are presented in Fig. [71 For the bulk modulus K both 
methods agree very well with the experiment and the 
difference between GGA and SE is small. For the shear 
modulus G, however, there is a significant difference be¬ 
tween GGA and SE. Fig. [8] shows the three sound veloc¬ 
ities vp, vb and vs- Both GGA and SE agree with the 
experiment pretty well, with vs (which is proportional 
to VG) displaying the largest SE-GGA difference. We 
also plot yj C44 /p (calculated with GGA+ZPV) as the 
golden dash-dot-dot curve in the lower panel of Fig. [51 
It has a meaning of vs calculated by using the modulus 
G44 instead of the averaged shear modulus G, like we did 
previously in Ref. E3- For an isotropic solid G = G44. 
For helium, using G44 instead of G underestimates vs 
by a few percent and worsens the agreement between 
GGA and the experiment. In other words, by calculat- 
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FIG. 8. (Color online) Aggregate sound velocities of hep he¬ 
lium versus pressure. 


ing vs from the proper Voigt-Reuss-Hill G in the present 
work we actually improve the GGA-experiment agree¬ 
ment compared to Ref. Hzl 

The bulk and shear moduli and sound velocities at TPa 
pressures are presented in Fig. [9j All these quantities 
show features at the metallization point, with kinks in G 
and, respectively, vs, being the most pronounced. 

The Debye temperature Tjj is presented in Fig. 0 
lower panel. Just like with shear modulus G, there is 
a significant difference between SE and GGA, with ex¬ 
perimental data points lying in the middle. The Debye 
temperature for the TPa pressures is plotted in Fig. [lOj 
upper panel. It shows a noticeable feature at the metal¬ 
lization point. 


D. Poisson’s ratio 

Figure [TT] shows the Poisson’s ratio cr as a function 
of pressure. The GGA+ZPV results^, obtained using 
G 44 instead of G, as explained above, are also plotted 
as the golden dash-dot-dot curve. In addition to the ex¬ 
perimental data of Zha. et alA, we also plot the low-p, 
low-T data of Nieto et al. 26 . Since the Poisson’s ratio 
([38]) is proportional to the difference 3 K — 2 G, it is a 
rather method-sensitive quantity. Just like for the elas¬ 
tic anisotropy parameters above, SE gives much stronger 
pressure dependence of a compared to GGA, and for the 
reasons outlined above for AP, ASi^ we find GGA re¬ 
sults more trustworthy in this aspect. On the other hand, 
the experimental pressure dependence of a seems to be 
closer to the SE one. Neither method agrees perfectly 
with the T = 300 K experimental data of Zha. et alA, 
however, both our method agree quantitatively with the 


FIG. 9. (Color online) Aggregate properties of hep helium in 
the TPa pressure range: sound velocities (upper panel), bulk 
and shear moduli (lower panel). Symbols are the calculated 
data points. The vertical line indicates the GGA metallization 
point. 



FIG. 10. (Color online) Debye temperature (upper panel) and 
Poisson’s ratio (lower panel) of hep helium in the terapascal 
pressure range. Symbols are the calculated data points. The 
vertical line indicates the GGA metallization point. 


p « 0, T « 0 result of Nieto et al— (a = 0.38, the purple 
dot in Fig. as long as ZPV is included. 

Although different theoretical and experimental meth¬ 
ods give somewhat different values of cr, they all agree 
upon one fact: da/dp < 0 , i.e. a decreases monotonously 
when pressure increases, fully confirming the surprising 
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result of Zha. et al In particular, while the role of ZPV 
for a is somewhat larger than for other quantities stud¬ 
ied in the present paper, and a shows an isotopic effect^, 
the negative pressure dependence of a is definitely not a 
quantum effect, as the ”No ZPV” curves show the same 
sign and order of magnitude of da / dp. This behavior of 
a is highly unusual, although solid hydrogen&I also has 
negative da/dp at least for low pressures. All heavier 
rare-gas solids (RGSs)£ have positive da/dp. The dif¬ 
ference in behavior between He and heavier RGSs is, we 
repeat, is not a quantum effect, i.e. it is not caused by the 
small mass of He atoms. Presumably the crucial factor 
here is the difference of the outermost electron shells in 
He (Is 2 shell) and heavier RGSs (ns 2 np 6 ). The situation 
is similar for c/a — ^/8/3, which have different sign for He 
and other RGS’s in GGA^ . In the language of SE poten¬ 
tial, both types of behavior of a can be accounted for by 
using pair and triple forces of exactly the same functional 
formas, w fth different values of parameters. 

The Poisson’s ratio at terapascal pressures is plotted 
in Fig. [lOj lower panel, a decreases monotonously up to 
highest pressures apart from a large feature at the met¬ 
allization point. There is no minimum in a (disregard¬ 
ing the step at the metallization point) and definitely no 
a —> 0.5 high-pressure asymptotic that many materials 
have. At highest pressures considered a is of the order of 
0.29. 

A question often asked is whether helium becomes a 
classical crystal under high pressure, or in other words, 
whether the relative effect of ZPV on various physi¬ 
cal quantities decreases with pressure. This question 
is surprisingly nontrivial, as it depends on the compli¬ 
cated interplay of the kinetic and potential energy of the 
He atoms. A recent diffusion Monte Carlo analysis^ 
has shown that the kinetic to potential energy ratio 
| Ekin/Epot | (which can be viewed as the measure of quan¬ 
tumness) does indeed decrease monotonuously with pres¬ 
sure, however this decrease becomes very slow for pres¬ 
sures p > 85 GPa. We have reached similar conclusions 
in the present work, which is best seen for Poisson’s ratio 
(fig. E]). a is dimensionless, and the difference between 
classical and quantum a slowly decreases with pressure. 
In particular, this difference is of the order of 0.01 in 
the 150 GPa pressure range, but in the TPa range (not 
shown) it reaches the value of about 0.0035. 

IV. CONCLUSION 

We have calculated five elastic constants of hep helium 
under pressure, and various derived quantities measured 
by Zha. et alA\ anisotropy parameters, sound velocities, 
Poisson’s ratio etc. We have analyzed these quantities 
both in the pressure range up to 150 GPa, where experi¬ 
mental data is available and semi-empirical potentials are 


applicable; and in the TPa pressure range (GGA only), 
where the metallization transition takes place. Both 



FIG. 11. (Color online) Poisson’s ratio of hep helium versus 
pressure. 


the experiment. Most calculated quantities display no¬ 
ticeable features at the metallization point. Zero-point 
vibrations do not affect elastic properties of helium in any 
dramatic way for the pressures considered in the present 
paper (disregarding the quantum crystal region at very 
low pressures). 

Our calculations predict significant elastic anisotropy 
for hep helium (A P ~ 1.14, A Si ~ 1.7, AS2 ~ 
0.93 at low pressures). Both GGA and SE overesti¬ 
mate anisotropy parameters compared to the experiment, 
which might be a temperature effect (experiments were 
carried out at T = 300 K). The three anisotropy param¬ 
eters become more isotropic (close to one) under TPa 
pressures, with the anisotropy of the shear modes Si and 
S 2 changing sign at the metallization point. Our calcu¬ 
lated Poisson’s ratio (PR) is in excellent agreement with 
the T = 0, p = 0 result of Ref. 0 (a = 0.38). Our 
calculations agree with the experimentally observed^ de¬ 
crease of the PR with pressure. Under TPa pressures PR 
reaches values ~ 0.31 at the metallization point (p ~ 17.5 
TPa) and ~ 0.29 at p = 30 TPa. We have shown that 
the negative sign of the pressure dependence of PR is 
not a quantum effect by performing calculations without 
zero-point vibrations, which yield similar results. 
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